Introduction

We wish to create conservation-sensitive, climate-smart fisheries closures for the ABNJ of the Pacific Ocean using prioritizr.

Prioritizr (Hanson et al., 2018) is a spatial prioritization software (R package) that uses integer linear programming (ILP) allowing it to create exact solutions to large spatial planning problems faster than its heuristic counterparts (e.g. Marxan).

This file runs multiple scripts found in the main Github repository: tinbuenafe/PacificProj, where the codes are explained in greater detail.

Assembling the data layers

Before running prioritizr, we first had to assemble the following data layers:

  1. Study Area
  2. Fisheries conservation features
  3. Climate smart features
  4. Cost layer

Study Area

The study area has the following characteristics:

  • Pacific ABNJ
  • Limited boundaries to areal jurisdiction of RFMOs responsible for managing large pelagic fisheries (IATTC and WCPFC)
  • Robinson projection
  • Divided into Longhurst provinces

We ran the following code using polygons from rnaturalearth package and the EEZ .shp files from Marine Regions (Flanders Marine Institute, 2019) to create a Pacific-centered map of the global ABNJs.

The relevant scripts for creating the study area are labeled with a prefix of 01 in ~/scripts/.

# to create the Pacific-centered maps:
source("scripts/01a_StudyArea_CreatingPUs.R")
# this script uses the following functions (embedded actual script):

# function to create polygon of the ABNJ by masking the EEZs and the landmasses
source("scripts/study_area/fCreateMaskedPolygon.R") 
# function to reproject these polygon to Robinson projection
source("scripts/study_area/fConvert2PacificRobinson.R") 
# function to create the hexagonal planning units
source("scripts/study_area/fCreate_PlanningUnits.R")

This figure shows the Pacific-centered global map (using Robinson projection).

world_robinson <- readRDS("outputs/01_StudyArea/01a_StudyArea/PacificCenterLand.rds")
global_plot <- ggplot() +
  geom_sf(data = world_robinson, colour = "grey20", fill="grey20", size = 0.1, show.legend = "line") +
  theme_bw()
global_plot

We then created equal-area hexagonal planning units, or PUs, (0.5\(^\circ\) x 0.5\(^\circ\) at the equator) of the study area. This resulted in 31,917 PUs.

# define objects for plotting; objects are the result of running 01_StudyArea
pacific_robinson <- readRDS("outputs/01_StudyArea/01a_StudyArea/PacificCenterABNJ.rds")
PUsPac <- readRDS("outputs/01_StudyArea/01a_StudyArea/PacificABNJGrid_05deg.rds")
Bndry <- fCreateRobinsonBoundary(west = 78, east = 140, north = 51, south = 60) # coordinates in degrees

study_area <- ggplot() +
    geom_sf(data = pacific_robinson, colour = NA, fill = NA, size = 0.2, show.legend = "line") +
    geom_sf(data = world_robinson, color = "grey20", fill="grey20", size = 0.1, show.legend = "line") +
    geom_sf(data = PUsPac, colour = "grey64", aes(fill = "ABNJ"), size = 0.1, show.legend = TRUE) + 
    scale_fill_manual(name = "Study Area", values = c("ABNJ" = "coral3")) +
    coord_sf(xlim = c(st_bbox(Bndry)$xmin, st_bbox(Bndry)$xmax), # Set limits based on Bndry bbox
            ylim = c(st_bbox(Bndry)$ymin, st_bbox(Bndry)$ymax),
            expand = TRUE) +
    theme_bw()
study_area

After demarcating the study area, we categorized the PUs into their Longhurst Provinces, .shp files from Marine Regions (Flanders Marine Institute, 2009), using the function fCategProv() found in this source code.

# function to assign the provinces defined in:
source("scripts/01b_StudyArea_fCategProv.R")
# running this function:
source("scripts/01c_StudyArea_fCategProvRun.R")
longhurst

Fisheries Conservation Features

The following were used in this study:

  1. Distribution of selected bycatch species
  2. Distribution of larval spawning grounds of selected commercially targeted species

Bycatch Distribution

First, we compiled a list of bycatch reported in the Pacific ABNJ. Then, we extracted their AquaMaps distribution. AquaMaps (Kaschner et al., 2019) is a niche-based model that predicts the relative occurrences of species from their environmental envelopes or their environmental and oceanographic preferences (Kesner-Reyes et al., 2016).

We limited the bycatch included in this study to sea turtles, but this workflow can be easily replicated and expanded to include more species.

The following were the sea turtles reported as bycatch and intersected with the study area:

  • Loggerhead turtle (Caretta caretta)
  • Green turtle (Chelonia mydas)
  • Leatherback turtle (Dermochelys coriacea)
  • Hawksbill turtle (Eretmochelys imbricata)
  • Olive ridley turtle (Lepidochelys olivacea)

The relevant scripts for creating the bycatch distribution layer are labeled with a prefix of 02_RawAQM and 03_AQM in ~/scripts/.

The plots below show the distribution of these species using a probability threshold of 0.5 and the functions fAquaStart() and fBycatchDistMap().

# to extract the distribution of all species that are included in the AquaMaps dataset:
# the function is found in:
source("scripts/02a_RawAQM_fAquaStart.R")
# the function is ran in:
source("scripts/02b_RawAWM_fAquaStartRun.R")

# to create the distribution maps of the sea turtle species included in this study:
# the function is found in:
source("scripts/02c_RawAQM_fBycatchDistMap.R")
# the function is ran in:
source("scripts/02d_RawAQM_fBycatchDistMapRun.R")

Spawning grounds of Commercially Targeted Species

Just like the bycatch layer, we compiled a list of commercially targeted large pelagic species managed in the Pacific from different sources. From this list, we have 4 species with available data on their global larval spawning grounds:

  • Yellowfin Tuna (Thunnus albacares)
  • Albacore (Thunnus alalunga)
  • Skipjack Tuna (Katsuwonus pelamis)
  • Swordfish (Xiphias gladius)

The models for these species are from a report by Mercer (2019) where he created Generalized Additive Models (GAMs) to predict spawning grounds. He used average seasonal abundance from a historical dataset of longline fisheries surveys (dated back to the 1950s), and oceanographic and environmental conditions from global databases.

We reran his models (see scripts labeled with a prefix 05a in ~/scripts/ for complete reruns) and used the medians of the predictions for each species as probability thresholds to create the distribution maps and the function fCommercialFeat().

# to intersect the results of the GAMs with the PUs
# the function is in:
source("scripts/05b_Commercial_fCommercialFeat.R")
# the function is ran at:
source("scripts/05c_Commercial_fCommercialFeatRun.R")
global_plots

Joining the conservation features’ layers with Longhurst provinces’ layer

Each province is expected to be biogeographically and ecologically unique. Including this layer in the spatial prioritization spreads the risk of exposure of populations to calamities, diseases or pollution. We joined the features’ layers with the provinces’ layer, creating multifaceted features, with each feature subtyped by provinces. This increased the number of features from 9 (5 from bycatch and 4 from large pelagics) to 66 features.

We used the function fProvIntersect().

# to intersect the conservation features with the provinces
source("scripts/03_Features_fProvIntersect.R")

# function ran with the bycatch layer:
source("scripts/03b_AQM_fProvIntersectRun.R")
# function ran with the commercially targeted species layer:
source("scripts/05e_Commercial_fProvIntersect.R")

This is an example of how the data looks like (so far) for bycatch:

head(bycatch)
## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3779074 ymin: 168708.9 xmax: -3695824 ymax: 473118.3
## CRS:           +proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
## # A tibble: 6 x 5
## # Groups:   province [1]
##   province cellsID species                                 geometry new_feature 
##   <chr>      <int> <chr>                              <POLYGON [m]> <chr>       
## 1 WARM           1 ITS-Mam-… ((-3751324 216773.5, -3779074 232795.… WARM_ITS-Ma…
## 2 WARM           2 ITS-Mam-… ((-3751324 312902.8, -3779074 328924.… WARM_ITS-Ma…
## 3 WARM           3 ITS-Mam-… ((-3751324 409032.1, -3779074 425053.… WARM_ITS-Ma…
## 4 WARM           4 ITS-Mam-… ((-3723574 168708.9, -3751324 184730.… WARM_ITS-Ma…
## 5 WARM           5 ITS-Mam-… ((-3723574 264838.2, -3751324 280859.… WARM_ITS-Ma…
## 6 WARM           6 ITS-Mam-… ((-3723574 360967.5, -3751324 376989,… WARM_ITS-Ma…

And for the large pelagics:

head(pelagic)
## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3779074 ymin: 216773.5 xmax: -3695824 ymax: 1482476
## CRS:           +proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
## # A tibble: 6 x 5
## # Groups:   province [2]
##   province cellsID species                                  geometry new_feature
##   <chr>      <int> <chr>                               <POLYGON [m]> <chr>      
## 1 WARM           1 ALB     ((-3751324 216773.5, -3779074 232795.1, … WARM_ALB   
## 2 WARM           2 ALB     ((-3751324 312902.8, -3779074 328924.4, … WARM_ALB   
## 3 WARM           3 ALB     ((-3751324 409032.1, -3779074 425053.6, … WARM_ALB   
## 4 WARM           5 ALB     ((-3723574 264838.2, -3751324 280859.7, … WARM_ALB   
## 5 WARM           6 ALB     ((-3723574 360967.5, -3751324 376989, -3… WARM_ALB   
## 6 NPSW           7 ALB     ((-3723574 1418389, -3751324 1434411, -3… NPSW_ALB

Each sf object has the following columns:

  • province is the longhurst province code
  • cellsID is the unique ID for each hexagonal PU
  • species is the species code (e.g. ALB = Albacore)
  • new_feature is the subtyped features code (e.g. WARM-ALB = Albacore in the Western Pacific Warm Pool)
  • geometry

Climate-smart Features

To create climate-smart closures, we prioritized protection of areas with:

  1. low levels of exposure to climate warming, and
  2. high levels of potential retention of biodiversity.

We determined these areas using the following metrics:

  1. Relative Climate Exposure (RCE) index

    \(\large RCE(yr^{-1}) = \frac{{\text Slope \hspace{2mm}} (^{\circ}Cyr^{-1} \hspace{2mm} 2050 - 2100)}{{\text Seasonal \hspace{2mm} Range} \hspace{2mm} (^{\circ}C \hspace{2mm} 2015 - 2020)}\)

  • numerator is the slope of the linear regression of projected mean annual temperatures
  • denominator is the current mean seasonal temperature range

Areas with low RCE index values are areas that will be exposed to low levels of warming relative to their current seasonal fluctuations of temperature. Areas of low RCE index values are prioritized for selection.

  1. Climate velocity

    \(\large Climate \hspace{2mm} velocity \hspace{2mm} (kmyr^{-1}) = \frac{{\text Slope \hspace{2mm}} (^{\circ}Cyr^{-1} \hspace{2mm} 2050 - 2100)}{{\text Spatial \hspace{2mm} gradient} \hspace{2mm} (^{\circ}Ckm^{-1} \hspace{2mm} 2050 - 2100)}\)

  • numerator is the slope of the linear regression of projected mean annual temperatures
  • denominator is calculated from the vector sum of the latitudinal and longitudinal pairwise differences of the mean temperature at each focal cell using a 3 x 3 neighborhood window

Areas of slow climate velocity are areas that are more likely to retain their current environmental conditions and consequentially, their biodiversity. Areas of slow climate velocity are prioritized for selection.

The temperatures used to calculate the aforementioned metrics were derived from 12 Global Circulation Models (GCMs) from the Couple Model Intercomparison Project Phase 6 (CMIP6; https://esgf-node.llnl.gov). The GCMs are under these three climate scenarios:

  1. SSP1-2.6: optimistic scenario
  2. SSP2-4.5: intermediate scenario
  3. SSP5-8.5: pessimistic scenario

RCE

We intersected the RCE values with the PUs using fClimateInt(). The relevant scripts for creating the climate-smart features layer are found in ~/scripts/ with the prefix 07.

# function is found at:
source("scripts/07a_Climate_fClimateInt.R")
# function is ran at:
source("scripts/07b_Climate_fClimateIntRun.R")

Climate velocity

We intersected the climate velocity values with the PUs using the same function, fClimateInt().

Cost layer

The cost layer represents the opportunity cost to fishing of establishing the proposed fisheries closures.

The cost layer for each planning unit was represented in US Dollars, and was calculated by multiplying the:

  • estimated total catch of medium - large pelagics in kg (from Watson, 2017),and
  • the ex-vessel price of the same species in the aforementioned functional groups in US$ \(kg^{-1}\) (from Tai et al., 2017)
# function is found at:
source("scripts/06a_Cost_fCostPU.R")
# function is ran at:
source("scripts/06b_Cost_fCostPURun.R")
## [1] "minimum: 0.0036497397813946 maximum: 5977.41455078125"

We used log-transformed cost for the figure below:

Representation Targets

Filtering and Retaining PUs

We join the climate-smart layer with the conservation features layer. Then, we only retain planning units that have RCE and/or climate-velocity values belonging to their lower 25th percentile range. This was how we priortized protection of areas that are predicted to have lower exposure to climate warming and higher retention of biodiversity.

We used the function fFilterQuartile() and the scripts that are relevant for this step are found in ~/scripts/ with the prefix 08.

# function is found at:
source("scripts/08a_Filter_fFilterQuartile.R")
# function is ran at:
source("scripts/08b_Filter_fFilterQuartile25Run.R")

This was done for both the bycatch and large pelagic layers, across different scenarios (SSP1-2.6, SSP2-4.5, SSP5-8.5). Below is an example of how the data looks like, with the following columns:

  • province: Longhurst provinces
  • cellsID: unique ID for each PU
  • species: species code
  • new_features: subtyped features (feature x province)
  • total_area: the total area covered by each subtyped feature (in km\(^{-1}\))
  • velocity: climate velocity
  • velo_tvalue: transformed values of climate velocity (velocity * 10)
  • RCE: RCE index values
  • RCE_tvalue: transformed values of RCE (cube root of RCE)
# for bycatch under the climate scenario SSP1-2.6.
bycatchSSP126 <- readRDS("outputs/08_Filter/08b_Filter25/bycatchSSP126_25percentile.rds") %>% 
  as_tibble() %>% 
  dplyr::select(-prov_descr, -area_km2, -cellsID.2, -geometry)
head(bycatchSSP126)
## # A tibble: 6 x 9
##   province cellsID species  new_features  total_area velocity velo_tvalue    RCE
##   <chr>      <int> <chr>    <chr>              <dbl>    <dbl>       <dbl>  <dbl>
## 1 WARM           1 Rep-343… WARM_Rep-343…   1664582.    2.16        21.6  0.0740
## 2 WARM           2 Rep-343… WARM_Rep-343…   1664582.    1.37        13.7  0.0626
## 3 WARM           3 Rep-343… WARM_Rep-343…   1664582.    0.531        5.31 0.0613
## 4 WARM           4 Rep-343… WARM_Rep-343…   1664582.    1.33        13.3  0.0765
## 5 WARM           5 Rep-343… WARM_Rep-343…   1664582.    2.91        29.1  0.0689
## 6 WARM           6 Rep-343… WARM_Rep-343…   1664582.    0.743        7.43 0.0585
## # … with 1 more variable: RCE_tvalue <dbl>
# for large pelagics under the cliamte scenario SSP5-8.5.
commercialSSP585 <- readRDS("outputs/08_Filter/08b_Filter25/commercialSSP585_25percentile.rds") %>% 
  as_tibble() %>% 
  dplyr::select(-prov_descr, -area_km2, -cellsID.2, -geometry)
head(commercialSSP585)
## # A tibble: 6 x 9
##   province cellsID species new_features total_area velocity velo_tvalue   RCE
##   <chr>      <int> <chr>   <chr>             <dbl>    <dbl>       <dbl> <dbl>
## 1 WARM           1 ALB     WARM_ALB        963004.     22.0        220.  3.47
## 2 WARM           3 ALB     WARM_ALB        963004.     25.7        257.  2.85
## 3 WARM           5 ALB     WARM_ALB        963004.     31.2        312.  3.18
## 4 WARM          14 ALB     WARM_ALB        963004.     23.4        234.  3.36
## 5 WARM          36 ALB     WARM_ALB        963004.     24.8        248.  3.25
## 6 WARM          56 ALB     WARM_ALB        963004.     26.2        262.  3.16
## # … with 1 more variable: RCE_tvalue <dbl>

Aside from these climate-smart closures, we also explored climate-uninformed, wherein the climate-smart features were not included in the analyses. Subsequently, the step of filtering and retention of planning units was also not included.

source("scripts/08c_Filter_fFilterQuartile100Run.R")

This is how the data looks like:

# for bycatch:
bycatch_uninformed <- readRDS("outputs/08_Filter/08c_Filter100/bycatchuninformed_100percentile.rds") %>% 
  dplyr::select(-prov_descr, -area_km2, -geometry) %>% 
  as_tibble()
head(bycatch_uninformed)
## # A tibble: 6 x 6
##   province cellsID species  new_features total_area                     geometry
##   <chr>      <int> <chr>    <chr>             <dbl>                <POLYGON [m]>
## 1 WARM           1 ITS-Mam… WARM_ITS-Ma…    912319. ((-3751324 216773.5, -37790…
## 2 WARM           2 ITS-Mam… WARM_ITS-Ma…    912319. ((-3751324 312902.8, -37790…
## 3 WARM           3 ITS-Mam… WARM_ITS-Ma…    912319. ((-3751324 409032.1, -37790…
## 4 WARM           4 ITS-Mam… WARM_ITS-Ma…    912319. ((-3723574 168708.9, -37513…
## 5 WARM           5 ITS-Mam… WARM_ITS-Ma…    912319. ((-3723574 264838.2, -37513…
## 6 WARM           6 ITS-Mam… WARM_ITS-Ma…    912319. ((-3723574 360967.5, -37513…
# for large pelagics:
commercial_uninformed <- readRDS("outputs/08_Filter/08c_Filter100/commercialuninformed_100percentile.rds") %>%
  dplyr::select(-prov_descr, -area_km2, -geometry) %>% 
  as_tibble()
head(commercial_uninformed)
## # A tibble: 6 x 6
##   province cellsID species new_features total_area                      geometry
##   <chr>      <int> <chr>   <chr>             <dbl>                 <POLYGON [m]>
## 1 WARM           1 ALB     WARM_ALB        963004. ((-3751324 216773.5, -377907…
## 2 WARM           2 ALB     WARM_ALB        963004. ((-3751324 312902.8, -377907…
## 3 WARM           3 ALB     WARM_ALB        963004. ((-3751324 409032.1, -377907…
## 4 WARM           5 ALB     WARM_ALB        963004. ((-3723574 264838.2, -375132…
## 5 WARM           6 ALB     WARM_ALB        963004. ((-3723574 360967.5, -375132…
## 6 NPSW           7 ALB     NPSW_ALB       4390870. ((-3723574 1418389, -3751324…

Assigning targets

To determine the relative targets per subtyped feature, we checked each one’s IUCN category (IUCN, 2015) using rredlist(). Threatened features (i.e. classified as vulnerable, endangered, threatened) were assigned a fixed maximum target. Features (\(i\)) that were not classified as threatened (e.g. classified as near-threatened, etc.) were assigned a relative target based on their distribution.

Relative target\(_i\) \(= Target_{max} \cdot \frac{PU_i}{PU_{total}} \cdot (Target_{max} - Target_{min})\)

  • PU_i: planning units of feature \(i\)
  • PU_{total}: total planning units
  • Target_{max}: Maximum target / the fixed maximum target assigned to threatened features (e.g. 1)
  • Target_{min}: Minimum target (e.g. 0.1)

Lower targets are assigned to features with broader distribution, and higher targets to features with more restricted distribution.

We reran this step for an array of targets, for both climate-smart and climate-uninformed closures:

Runs Maximum Minimum
1 1 0.1
2 0.9 0.1
3 0.8 0.1
4 0.7 0.1
5 0.6 0.1
6 0.5 0.1
7 0.4 0.1
8 0.3 0.1
9 0.2 0.1
10 0.1 0

We used the function fRepresentTarget(), and the scripts are found in ~/scripts/ with the prefix 09.

# the function is found in:
source("scripts/09a_Target_fRepresentTarget.R")
# the runs for climate-smart closures are found in:
source("scripts/09b_Target_fRepresentTargetSmartRun.R")
# the runs for climate-uninformed closures are found in:
source("scripts/09c_Target_fRepresentTargetUninformedRun.R")

For this .rmd, we’ll be showing results using Runs 2 (Max = 0.9, Min = 0.1), 5 (Max = 0.6, Min = 0.1) and 8 (Max = 0.3, Min = 0.1). This workflow shows that it can be replicated using other parameters.

Below is an example of how the data looks like:

target90_commercial <- readRDS("outputs/09_Target/09b-c_TargetRuns/02_Target90/SSP126/target_commercialSSP126.rds") %>% 
  as_tibble()
head(target90_commercial)
## # A tibble: 6 x 5
##   new_features    total_area category                            geometry target
##   <chr>                <dbl> <chr>                         <GEOMETRY [m]>  <dbl>
## 1 ARCH_ALB           218743. NT       POLYGON ((-2086316 -2859363, -2114… 0.185 
## 2 ARCH_SKP            18673. LC       MULTIPOLYGON (((-2086316 -2859363,… 0.0158
## 3 ARCH_SWO            32011. LC       MULTIPOLYGON (((-2086316 -2859363,… 0.0271
## 4 CHIL_YFT            29344. NT       MULTIPOLYGON (((8985986 -504196, 8… 0.0248
## 5 non-categ_Long…    424148. NT       MULTIPOLYGON (((-476808.4 -2859363… 0.359 
## 6 non-categ_Long…    242752. LC       MULTIPOLYGON (((-698809.5 -1705812… 0.205
target90_bycatch <- readRDS("outputs/09_Target/09b-c_TargetRuns/02_Target90/SSP126/target_bycatchSSP126.rds") %>% 
  as_tibble()
head(target90_bycatch)
## # A tibble: 6 x 5
##   new_features   total_area category                             geometry target
##   <chr>               <dbl> <chr>                          <GEOMETRY [m]>  <dbl>
## 1 ARCH_Rep-1347…     88031. VU       POLYGON ((-1753314 -2859363, -17810…     90
## 2 ARCH_Rep-2226…     88031. EN       POLYGON ((-1753314 -2859363, -17810…     90
## 3 ARCH_Rep-3437…    445489. VU       MULTIPOLYGON (((-726559.6 -3003557,…     90
## 4 ARCH_Rep-4127…     96034. CR       POLYGON ((-1753314 -2859363, -17810…     90
## 5 ARCH_Rep-5414…     13338. VU       POLYGON ((-1864315 -2763234, -18920…     90
## 6 CCAL_Rep-3437…    554861. VU       MULTIPOLYGON (((5794721 2139359, 57…     90
  • new_features: subtyped features (feature x province)
  • total_area: the total area covered by each subtyped feature (in km\(^{-1}\))
  • category: IUCN redlist category
  • target: assigned target in percentage

Running prioritizr

prioritizr (v. 7.0.1.2) uses the following functions to define the problem.

# installing latest version of prioritizr from GitHub
devtools::install_github("prioritizr/prioritizr")
# installing gurobi as the optimizer
install.packages("/Library/gurobi911/mac64/R/gurobi_9.1-1_R_4.0.2.tgz", repos = NULL)
library(priortizr)
library(gurobi)

p1 <- prioritizr::problem(x, feature_list, "cost") %>% 
    add_min_set_objective() %>%
    add_relative_targets(target_list) %>%
    add_binary_decisions() %>%
    add_gurobi_solver(gap = 0.1, verbose = FALSE) 

s1 <- solve(p1)
  • problem(): we define the data
  • x: data frame having all the data layers (all planning units with their specific ID (cellsID), features data, cost)
  • feature_list: list of the features in x
  • cost: name of the column in x where the cost is found
  • add_min_set_objective(): function that allows us to employ the minimum set objective function. This allows prioritizr to create spatial plans that will minimize the cost.
  • add_relative_targets(): target_list list of assigned targets for each feature
  • add_binary_decisions(): prioritizr will only either select or not select a PU
  • add_gurobi_solver(): Gurobi was used as the optimizer to solve the problem (at a default 10%)
  • solve(): prioritizr solves the problem p1

Prioritizr Runs

We used the fPrioritizrRun() function to streamline the runs. fPrioritizrRun() arranges all the data and makes sure that they’re in the right format, then runs problem() and solve(). This function also saves the summaries each plan (cost, % area protected and number of representative features that achieved 30% or more protection).

The relevant scripts for this step is found in ~/scripts/ with the prefix 10.

# the function is in:
source("scrpits/10a_Prioritizr_fPrioritizr.R")
# to run the function with climate-smart data
source("scripts/10b_Prioritizr_fPrioritizrAQMSmartRun.R")
# to run the function with climate-uninformed data
source("scripts/10d_Prioritizr_fPrioritizrAQMUninformedRun.R")
# to create no-regret closures with the climate-smart plans

To create no-regret plans, we used thre function fCreateNoRegret() which intersected the selected PUs in all three climate scenarios. The relevant scripts for this step is found in ~/scripts with the prefix 11.

# the function is in:
source("scripts/11a_NoRegret_fCreateNoRegret.R")
# the runs are found in:
source("scripts/11b_NoRegret_fCreateNoRegretAQMRuns.R")

As mentioned in the previous section, we will only be showing results of Runs 2 (Max = 0.9, Min = 0.1), 5 (Max = 0.6, Min = 0.1) and 8 (Max = 0.3, Min = 0.1).

Run 2: Maximum Target = 90%

Climate-smart runs

SSP1-2.6 run

SSP126_target90 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/02_Target90/solution_SSP126.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario    cost perc_area_protected protected_PU
##    <dbl> <chr>      <dbl>               <dbl>        <dbl>
## 1     90 SSP1.26  267943.                26.0           59

SSP2-4.5 run

SSP245_target90 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/02_Target90/solution_SSP245.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario    cost perc_area_protected protected_PU
##    <dbl> <chr>      <dbl>               <dbl>        <dbl>
## 1     90 SSP2.45  267309.                26.9           61

SSP5-8.5 run

SSP585_target90 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/02_Target90/solution_SSP585.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario    cost perc_area_protected protected_PU
##    <dbl> <chr>      <dbl>               <dbl>        <dbl>
## 1     90 SSP5.85  313651.                29.8           60

No-regret plan

noregret_target90 <- readRDS("outputs/11_NoRegret/11b_AQMRuns/noregretclosures_Target90.rds")

Climate-uninformed run

uninformed_target90 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/02_Target90/solution_uninformed.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario      cost perc_area_protected protected_PU
##    <dbl> <chr>        <dbl>               <dbl>        <dbl>
## 1     90 uninformed 744951.                73.5           66

Run 5: Maximum Target = 60%

Climate-smart runs

SSP1-2.6 run

SSP126_target60 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/05_Target60/solution_SSP126.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario    cost perc_area_protected protected_PU
##    <dbl> <chr>      <dbl>               <dbl>        <dbl>
## 1     60 SSP1.26  132371.                17.3           54

SSP2-4.5 run

SSP245_target60 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/05_Target60/solution_SSP245.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario    cost perc_area_protected protected_PU
##    <dbl> <chr>      <dbl>               <dbl>        <dbl>
## 1     60 SSP2.45  131116.                18.0           54

SSP5-8.5 run

SSP585_target60 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/05_Target60/solution_SSP585.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario    cost perc_area_protected protected_PU
##    <dbl> <chr>      <dbl>               <dbl>        <dbl>
## 1     60 SSP5.85  155441.                19.9           55

No-regret plan

noregret_target60 <- readRDS("outputs/11_NoRegret/11b_AQMRuns/noregretclosures_Target60.rds")

Climate-uninformed run

uninformed_target60 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/05_Target60/solution_uninformed.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario      cost perc_area_protected protected_PU
##    <dbl> <chr>        <dbl>               <dbl>        <dbl>
## 1     60 uninformed 369968.                49.0           63

Run 8: Maximum Target = 30%

Climate-smart runs

SSP1-2.6 run

SSP126_target30 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/08_Target30/solution_SSP126.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario   cost perc_area_protected protected_PU
##    <dbl> <chr>     <dbl>               <dbl>        <dbl>
## 1     30 SSP1.26  49676.                8.69           39

SSP2-4.5 run

SSP245_target30 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/08_Target30/solution_SSP245.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario   cost perc_area_protected protected_PU
##    <dbl> <chr>     <dbl>               <dbl>        <dbl>
## 1     30 SSP2.45  49934.                9.01           39

SSP5-8.5 run

SSP585_target30 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/08_Target30/solution_SSP585.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario   cost perc_area_protected protected_PU
##    <dbl> <chr>     <dbl>               <dbl>        <dbl>
## 1     30 SSP5.85  57173.                9.96           38

No-regret plan

noregret_target60 <- readRDS("outputs/11_NoRegret/11b_AQMRuns/noregretclosures_Target60.rds")

Climate-uninformed run

uninformed_target30 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/08_Target30/solution_uninformed.rds")

  • cost is in USD
  • perc_area_protected is the percentage of the Pacific ABNJ protected
  • protected_PU is the number of features that have >= 30% (by area) protection
## # A tibble: 1 x 5
##   target scenario      cost perc_area_protected protected_PU
##    <dbl> <chr>        <dbl>               <dbl>        <dbl>
## 1     30 uninformed 129125.                24.5           42

Summary

Since we ran prioritizr from 10 - 100% maximum target representation (in increments of 10%), we’re showing a summary of how cost and percent protection varied with increasing target representation.

source("scripts/12c_Stat_SomeAnalyses.R")
# log-transformed cost with increasing target representation:
log_plot

# % area protected with increasing target representation:
area_plot

# number of planning units that achieved 30% or more protection (maximum is 66)
protected_features

We also determined Cohen’s Kappa Coefficient between plans created for each run. Below are the summaries for each run as well as a matrix plot of the correlations between plans.

The relevant scripts for this part are found in ~/scripts/ with prefix 12.

# the function for creating the correlation plots is found in:
source("scripts/12d_Stat_fKappaCorrplot.R")
# the runs are found in:
source("scripts/12e_Stat_fKappaCorrplotRun.R")

Summary for Run 2: Maximum Target = 90%

## # A tibble: 4 x 5
##   target scenario      cost perc_area_protected protected_PU
##    <dbl> <chr>        <dbl>               <dbl>        <dbl>
## 1     90 SSP1.26    267943.                26.0           59
## 2     90 SSP2.45    267309.                26.9           61
## 3     90 SSP5.85    313651.                29.8           60
## 4     90 uninformed 744951.                73.5           66

Summary for Run 5: Maximum Target = 60%

## # A tibble: 4 x 5
##   target scenario      cost perc_area_protected protected_PU
##    <dbl> <chr>        <dbl>               <dbl>        <dbl>
## 1     60 SSP1.26    132371.                17.3           54
## 2     60 SSP2.45    131116.                18.0           54
## 3     60 SSP5.85    155441.                19.9           55
## 4     60 uninformed 369968.                49.0           63

Summary for Run 8: Maximum Target = 30%

## # A tibble: 4 x 5
##   target scenario      cost perc_area_protected protected_PU
##    <dbl> <chr>        <dbl>               <dbl>        <dbl>
## 1     30 SSP1.26     49676.                8.69           39
## 2     30 SSP2.45     49934.                9.01           39
## 3     30 SSP5.85     57173.                9.96           38
## 4     30 uninformed 129125.               24.5            42